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Abstract 

We simulate the dynamics of phase assembly in binary immiscible fluids and ternary microemulsions 
using a three-dimensional hydrodynamic lattice gas approach. For critical spinodal decomposition we 
perform the scaling analysis in reduced variables introduced by [Q, ^] . We find a late-stage scaling 
exponent consistent with the R ~ t'^''^ inertial regime. However, as observed with the previous lattice- 
gas model of Appert et al. ^ our data does not fall in the same range of reduced length and time as 
that of Kendon et al. QQ. For off-critical binary spinodal decomposition we observe a reduction of 
the effective exponent with the volume fraction of the minority phase. However, the n — ^ Lifshitz- 
Slyzov- Wagner droplet coalescence exponent is not observed. Adding a sufficient number of surfactant 
particles to a critical quench of binary immiscible fluids produces a ternary bicontinuous microemulsion. 
We observe a change in scaling behaviour from algebraic to logarithmic growth for amphiphilic fluids in 
which the domain growth is not arrested. For formation of a microemulsion where the domain growth is 
halted we find a stretched exponential growth law provides the best fit to the data. 

1 Introduction 

The study of phase ordering kinetics has become a testbed for new complex fluid simulation methods. Despite 
intensive analysis by many methods, it remains a field in which many interesting and fundamental questions 
go unanswered. Constructing a model which correctly includes hydrodynamics but which is computationally 
simple enough to reach late times is a major theoretical challenge. New mcsoscale models such as lattice gas 
(LG) lattice Boltzmann (LB) dissipative particle dynamics (DPD) and Boltzmann-Vlasov |^ 

treatments have been crucial in increasing our understanding of these systems over the past decade. Such 
methods require significantly smaller computational resources than earlier molecular dynamics (MD) and 
Cahn-Hilliard approaches, and can therefore access more easily the late time regime in which hydrodynamics 
plays an important role. 

A much studied system which exhibits hydrodynamic influences on phase segregation is a f : f mixture 
of immiscible fluids quenched into the two phase region of its phase diagram (a so-called critical quench). 
Such a system undergoes spinodal decomposition, where initially bicontinuous domain structure coarsens 
by surface-tension driven flows. Spinodal decomposition elicits much interest because of its fundamental 
and industrial importance. For example, the mechanical properties of alloys depend on the dynamics of 
the phase separation process. Despite extensive theoretical [0-|9|, numerical ||l|-|^,[lO|-p^ and experimental 
investigation []l9| , [20[| some doubt remains about the true asymptotic late time growth behaviour of these 
systems. 
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The dynamics of phase segregation in binary immiscible fluids are also dependent on the composition of 
the mixture. Mixtures which do not have a 1 : 1 ratio of the species (so-called off-critical quenches) arc much 
less studied than their critical counterparts. As the quantity of the minority phase decreases, the domain 
structure ceases to be bicontinuous, so that nucleation and coalescence of droplets dominate the coarsening 
mechanism. 

The addition of surfactant to a system of binary immiscible fluids radically alters the equilibrium prop- 
erties and growth dynamics of such mixtures ||2l|] . In particular, the adsorption of surfactant (amphiphilic) 
molecules at oil-water interfaces leads to a dramatic reduction in the intcrfacial tension. This property is 
the origin of much industrial interest in surfactant systems. Only recently have computational techniques 
and sufflciently powerful parallel platforms become available which permit numerical simulation of the hy- 
drodynamic behaviour of amphiphilic fluids in three dimensions. Our hydrodynamic lattice gas model has 
recently been implemented in three dimensions, and its ability to reproduce many important phenomenolog- 
ical features of amphiphilic systems confirmed p^ . In the present paper we demonstrate this model's ability 
to capture quantitative features of the dynamics of binary immiscible and ternary amphiphilic fluids. 

Continuum approaches regard spinodal decomposition as a solution of the equations of motion for two 
phase flow. These are given by the Navier-Stokes equation within each phase, and by boundary conditions 
at the interfaces. The incompressible Navier-Stokes equations are, for each fluid phase: 



V • V = (2) 

where v is the fluid velocity, 77 is the shear viscosity and p is the scalar pressure. The first boundary condition 
is that the velocity of the two phases must match that of the interface: 

ui • n = U2 ■ n = Mint ■ n (3) 

where n is the interface normal and Ui, Ui and Ui„t are respectively the velocities of phase 1 and 2 at the 
interface, and the velocity of the interface itself. The second boundary condition requires that the stress 
difference at the interface is balanced by the interfacial tension: 

T..„-T,^„.,(i- + i-)„ (4) 
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where T-^ and denote the stress tensor in phases 1 and 2, Ra and Rb are the two principal radii of the 
interfaces and a is the interfacial tension. 

Clearly, a numerical integration scheme for the above equations would be of enormous computational 
complexity. In order to obtain a numerical model which captures the essential physics in a computationally 
tractable way, an alternative approach is required. Remaining with a phenomenological description of the 
fluid system, it is convenient to define a scalar order parameter ip: 



V'(x) = Po(x) - /9^(x) 



(5) 



where /5o(x) and Pwi'x.) are the densities of oil and water at position x. One may then write a Cahn-Hilliard 
equation for the evolution of the order parameter: 



'dt 



■ V • VV' = AV> + AS,, 



(6) 



and. 



and 



F= d'^ 



SF 



(7) 



(8) 



This equation is coupled to the Navier-Stokes equations through a forcing term proportional to the gradient 
of the chemical potential: 

p + V • Vv^ = j/V^v - Vp - i/'V// + Av. (9) 

where and Av are gaussian noise fields satisfying an appropriate fluctuation-dissipation theorem . 

The equations - ^ define a model for immiscible fluid dynamics commonly referred to as Model H. 
This model enables one to obtain the scaling regimes for critical spinodal decomposition by dimensional 
analysis. The viscous regime is obtained when one may neglect the inertial terms on the left hand side of 
equation (m; then, balancing the forcing terms by the viscous terms one obtains: 



R{t) 1 



(10) 
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The inertial regime is obtained by balancing the forcing terms ipS/iJ by the inertial terms in eqn (|l|) : 

R(t) ^ (^)l/3i2/3 
P 

Equating these two regimes imphes that the crossover from viscous to inertial scaling occurs at i? ~ rf' /ap. 
However, no three-dimensional simulation method so far developed can reach a sufficient range of length and 
timescales to observe both viscous and inertial regimes in a single simulation. The work of Jury et al. Q 
and Kendon et al. overcomes this difficulty by introducing scaling variables ^-nd Tg which enable data 
from separate simulations to be combined. 

If the only physical effects involved in critical spinodal decomposition are capillary forces, viscous dis- 
sipation and fluid inertia, then the parameters governing domain growth are the surface tension cr, fluid 
mass density p, and viscosity rj. As emphasised by Jury et al. only one length, io ~ I P'^^ s-^d one 
time To — if' j pa^ can be constructed from these parameters. Data sets (L,T) from any model of spinodal 
decomposition can be expressed in units of reduced length: 

I^LIU, (12) 

and time, 

t = T/To. (13) 
If all other physics is excluded from late-stage growth, then the dynamical scaling hypothesis 1^ states that: 

lit)r.a + fit), (14) 



where l{t) is the domain size expressed in reduced length units and a is a non- universal constant which allows 
for a period of growth in which molecular diffusion leads to the formation of sharp interfaces. The function 
f{t) should then approach a universal form, identical for all 50 : 50 incompressible binary fluid mixtures 
following a deep quench. 

The Model H equations have been integrated numerically in three dimensions both with A ^ {a quench 
to finite temperature), and with A = (a quench to zero temperature) by a number of integration schemes 
The free energy defined in eqn (^) is also the basis for the lattice-Boltzmann method of Kendon et 
al. Q . Although the lattice-Boltzmann method has a physical motivation which replaces velocity and density 
fields by single particle distribution functions, it is essentially equivalent to direct numerical integration of 
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Model H such as that performed in ||12|,|l^] (i.e. it should be regarded as a finite-difference scheme for solving 
these equations in the absence of noise). The work of [^|l2|,|l^ is therefore a confirmation of the scaling 
laws derived above. The contribution of Kcndon et al. in whose work both growth regimes were clearly 
seen for the first time, shows that the simple arguments given above are incorrect in the crossover regime, 
and that the crossover in Model H extends over 10^ > t > 10® in reduced time imits. It should be noted 
that this crossover was previously believed from simple scaling arguments to occur at t = 1. 

Recently there has been further theoretical work by Grant and Elder: those authors suggest ||2^ that 
growth in the incrtial regime implies a Reynolds number that increases without limit. This would eventually 
lead to turbulent remixing of the fluids; the requirement that the Reynolds number be self-limiting (i.e. 
that the critical Reynolds number is not reached and therefore that turbulent remixing docs not occur) in 
the asymptotic regime implies a growth exponent of < 1/2. There is at present no numerical evidence for 
the n < 1/2, and recent theoretical challenges to this idea have been made. Novik and Coveney point out 
in that the relative strength of the interface (characterised by the Weber number) must also be taken 
into account. A small Weber number could delay the onset of turbulent remixing indefinitely. 

The hydrodynamic lattice gas used in the present paper does not rest on the macroscopic free energy 
functional proposed in cqn (^. The model, which is particulate in nature, is described in detail below, 
and has a theoretical justification from a 'bottom up' perspective. In a bottom-up view the lattice-gas 
technique may be regarded as a simplification of the molecular dynamics of a binary fiuid, abstracting the 
key microscopic properties in a fictitious microworld. Recent work has derived a microscopic basis for the 
Rothman-Kcller model of binary immiscible fiuids |^,|2^. However, no systematic method exists for coarse- 
graining a real molecular dynamics description of a system to the lattice-gas model used here (although such 
a systematic method does now exist for the DPD algorithm [2^]). 

In our model equations (0^) are emergent macroscopic properties. The single phase FHP lattice gas 
is known to satisfy eqn. (pi), and far from interfaces this behaviour is reproduced in our model l27LpO[. In 



Section 3.2 we demonstrate that our model has realistic surface tension behaviour at interfaces. 

The purpose of the present paper is to investigate the dynamics of domain growth of both binary im- 
miscible and ternary amphiphilic phases in our model. Section ^ contains a brief description of our model, 
while Section ^ contains a description of the quantitative measurements of surface tension and viscosity. In 
Sections ^, ^ and ^ we present our results for self-assembly in critical and off-critical binary immiscible and 
ternary amphiphilic fiuids respectively. We close the paper with discussion and conclusions in Section 0. 
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2 The lattice gas model 

Our lattice-gas model is based on a microscopic, bottom-up approach, where dipolar amphipile particles are 
included alongside the immiscible oil and water species. Lattice-gas particles can have velocities C;, where 
I < i < b, and b is the number of velocities per site. We shall measure discrete time in units of one lattice 
timcstep, so that a particle emerging from a collision at site x and time t with velocity Cj will advect to site 
X + Ci where it may undergo the next collision. We let nf (x, t) e {0, 1} denote the presence (1) or absence 
(0) of a particle of species a G {R,B,A} {R, B, A denoting red (oil), blue (water) and green (amphiphilc) 
species respectively) with velocity c.^, at lattice site x G £ and time step t. The collection of all nf (x, t) for 
1 < i < & will be called the population state of the site; it is denoted by 

n(x,i)GA/' (15) 

where we have introduced the notation JV for the (finite) set of all distinct population states. The amphiphile 
particles also have an orientation denoted by cr^ (x, t). This orientation vector, which has fixed magnitude a, 
specifies the orientation of the amphiphile particle at site x and time step t with velocity c^. The collection 
of the b vectors cr^ (x, t) at a given site x and time step t is called the orientation state. We also introduce 
the colour charge associated with a given site, 

g,(x,t) = nf (x,t)-7if (x,i), (16) 

as well as the total colour charge at a site, 

b 

9(x,t) = ^g,(x,0. (17) 

i=l 

The state of the model at site x and time t is completely specified by the population state and orientation 
state of all the sites. The time evolution of the system is an alternation between an advcction or propagation 
step and a collision step. In the first of these, the particles move in the direction of their velocity vectors to 
new lattice sites. This is described mathematically by the replacements 

<(x + c„t + l)^7if (x,i), (18) 

CTi (X + C,;,t + 1) ^ (T, (x,i) , (19) 
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for all X G £, I < i < b and a £ {R, B, A}. That is, particles with velocity simply move from point x 
to point X + Ci in one time step. In the collision step, the newly arrived particles interact, resulting in new 
momenta and surfactant orientations. The coUisional change in the state at a lattice site x is required to 
conserve the mass of each species present 

b 

p"(x,i)^5]<(x,t), (20) 

i 

as well as the £)-dimensional momentum vector 

6 

p(x,i)^^^c,<(x,t), (21) 

a i 

(where we have assumed for simplicity that the particles all carry unit mass). Thus, the set TV of population 
states at each site is partitioned into equivalence classes of population states having the same values of these 
conserved quantities. 

We assume that the characteristic time for collisional and orientational relaxation is sufficiently fast in 
comparison to that of the propagation that we can model this probability density as the Gibbsian equilibrium 
corresponding to a local, sitewise Hamiltonian function; that is 

Vis')^^cxp[^f3H{s')], (22) 

where /3 is an inverse temperature, H{s') is the energy associated with collision outcome s', and Z is the 
equivalence-class partition function. The sitewise Hamiltonian function for our model has been previously 
derived and described in detail for the two-dimensional version of the model |]3l[|, and we use the same one 
here. It is 

H{s') = J ■ (aE + ^P) + cr' ■ (eE + (P) + J : {e£ + (V) + -v(x, t)\ (23) 
where we have introduced the colour flux vector of an outgoing state, 

b 

J(x,t)=^c,qax,t), (24) 

1=1 

the total director of a site, 

b 

a (x, t)=Y, (x, i) ■ (25) 
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the dipolar flux tensor of an outgoing state, 



J(x,t) = ^c,a^ (x,t), 



(26) 



the colour field vector, 



the dipolar field vector. 



the colour field gradient tensor. 



the dipolar field gradient tensor, 



E(x,t) EE ^Cig(x + c,,i) , 



P(x,t) = -^c,S'(x + c,,t), 



f (x,t) =^c,E(x + c„i), 



V {x,t) = - '^ c^CiS (x + Cj, t) . 

i=l 

defined in terms of the scalar director field, 

b 

5(x,t)=^C, •^T,(x,t). 
i=l 

and the kinetic energy of the particles at a site. 



^|v(x,t)r, 



(27) 



(28) 



(29) 



(30) 



(31) 



(32) 



where v is the average velocity of all particles at a site, the mass of the particles is taken as unity, and a, ^, 
e, and S are coupling constants. To maintain consistency with previous work we use the coupling constants 
as previously defined in |^ . The values of these constants are: 

a = 1.0, e = 2.0, fi = 0.75, C = 0.5 

These values were chosen in order to maximise the desired behaviour of sending surfactant to oil-water 
interfaces while retaining the necessary miccllar binary water-surfactant phases. It should be noted that the 
inverse temperature-like parameter (3 (whose numerical value is varied in this paper) is not related in the 
conventional way to the kinetic energy. For a discussion of the introduction of this parameter into lattice 
gases we refer the reader to the original work by Chen, Chen, Doolen and Lee and Chan and Liang [^ . 
Eqs. (|2^)-(|3C|) were derived by assuming that there is an interaction potential between colour charges, and 
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that the surfactant particles are hke "colour dipoles" in this context pjl| . The term parameterised by a 
models the interaction of colour charges with surrounding colour charges as in the original Rothman-KcUer 
model that parameterised by ^ describes the interaction of colour charges with surrounding colour 
dipoles; that parameterised by e accounts for the interaction of colour dipoles with surrounding colour 
charges (alignment of surfactant molecules across oil- water interfaces); and finally that parameterised by 
C describes the interaction of colour dipoles with surrounding colour dipoles (corresponding to interfacial 
bending energy or "stiffness"). This model has been extensively studied in two dimensions ^4|-^, 
and the three-dimensional implementation employed in the present paper is described in more detail by 
Boghosian, Coveney and Love |22| . 

3 Viscosity and surface tension 

In this section we present the results of simulations designed to measure the values of the macroscopic physical 
parameters which control domain growth in fluid systems in our model, according to the top-down theories 
described in Section 1. As emphasised above, these parameters are the surface tension cr, the viscosity rj and 
the density p. 

3.1 Viscosity measurements 

We measured the viscosity by analysing the decay of shear waves. We performed simulations to observe the 
decay of shear waves with an initial velocity profile of the form: 



where is the lattice size in the x direction, vo{t) is the amplitude of the shear wave at time t, and is 
the imit vector in the z direction. Solving the Navier-Stokes equations for the time evolution of wo(i) gives: 



Where wo(0) is the initial amplitude of the shear wave and v is the kinematic viscosity. We therefore initialise 
the system with a velocity profile given by (|33| ) and observe the decay of the shear wave by calculating: 



v(x,t) = wo(0cos(27r— )e^ 



(33) 



X 



z;o(t)-«o(0)exp(-K^)2t). 



(34) 




(35) 
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where X^yz indicates summation over all lattice sites in the yz plane. This gives the mean z component of 
velocity, but includes any net z velocity the system may possess. We subtract this to obtain the velocity due 
only to the shear wave, 

y(x,i) = \4-i-^K(a^,t), (36) 

X 

and calculate V{){t) as the r.m.s. value of this quantity, 



vo{t)^ ^^Y.[V{x,t)Y. (37) 

We performed simulations at 5 different amplitudes of initial velocity profile and obtained kinematic viscosity 
V = 0.78 ± 0.01 in lattice units. 

3.2 Surface tension analysis 

In the present section we analyse the surface tension behaviour in our model for both binary immiscible 
and ternary amphiphilic systems. The central feature of ternary amphiphilic fluids is the lowering of the 
interfacial tension between oil and water by the adsorption of surfactant at the interface. The existence of 
a spinodal point in a binary immiscible fluid is an important feature of the fluid's thermal behaviour. At 
temperatures above the spinodal point the fluid will not demix into single phase domains. As one increases 
the temperature towards the spinodal point the surface tension should be reduced to zero. 

We used a direct dynamical method for calculating the surface tension across a flat interface. In the 
vicinity of an interface the pressure is locally anisotropic, as the pressure in the direction parallel to the 
interface is reduced by the tension on the interface itself. For a system with flat interface perpendicular to 
the z axis the surface tension is given by the line integral over z of the component Pn{z) of pressure normal 
to the interface minus the component Pt{z) tangential to the interface ||3^ : 

{PN{z)~PT{z))dz. (38) 

J —oo 

Where: 

Pt = Pxx = Pyy (39) 

Pn - Pzz (40) 

We begin by showing that our model is capable of reproducing physically realistic interfacial tensions in 
systems of binary immiscible fluids. We performed 12 simulations for values of 0.02 > /3 > 100, using systems 



11 



of size 64"^. The surface tension was measured every time step for 1000 time steps . The simulation with 
P — 0.02 was above the spinodal point. The spinodal point itself was located by a computational steering 
technique, which was found to be computationally more efficient than traditional taskfarm methods. The 
value of /3 was modified during a simulation and the phase separation behaviour visualised in real time. The 
value of (3 at the spinodal point was foimd to be 0.025 ± 0.003. This 'compusteering' technique represents a 
powerful and economic simulation technique, and will be the subject of a future paper. Equilibration effects 
were found to be significant only close to the spinodal point, and equilibration times for (3 ~ 0.03, 0.04, 0.05 
were taken as 400, 200 and 200 time steps respectively. All other simulations were allowed to equilibrate 
for 50 time steps and the surface tension was then time averaged for the remainder of the simulation. The 
surface tension as a function of (i is shown in fig. (|l]). 

We next analyse the behaviour of the surface tension as a function of surfactant concentration adsorbed 
onto an oil-water interface. We use two types of system, both initialised with a flat oil-water interface 
perpendicular to the z direction. The first system is initialised with a fixed amount of surfactant at the 
interface. The second is initialised with the surfactant dispersed throughout the system. Wc compute the 
surface tension as above from eqn. (|38|). 

The equilibrium state for such a system in our model is quite complex. The simplistic macroscopic 
picture of interfaces in ternary systems is fairly static, with a monolayer of surfactant coating the interface 
and lowering the surface tension. However, in reality as well as in our model the surfactant may also exist 
in bulk solution far from the interface, either as monomers or micelles or both. Equilibration of the system 
is achieved when the rates for adsorption and desorption of monomers and micelles are equal. These times 
are typically long, varying from 10000 to 20000 time steps in our model. 

For the first type of system we simulate, where all surfactant particles initially reside at the interface, 
the surface tension prior to equilibration calculated by ( |38| ) is an increasing function of time. The surfactant 
density at the interface decreases with time as monomers and micelles desorb into the bulk. For the second 
type of system, where all surfactant particles reside in the bulk, the surfactant density at the interface 
increases as the surfactant adsorbs to the interface. The surface tension calculated by (^8|) is a decreasing 
function of time. 

We calculate the surface tension as a function of surfactant density in two ways. For low surfactant 
concentrations where we can reach timescales on which the system is equilibrated we use the surface tension 
and intcrfacial density for the equilibrated system. For higher surfactant densities where the equilibration 
times become prohibitive wc plot the surface tension at time t against the density at time t. The results are 
shown in fig. (||). 
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4 Critical Spinodal Decomposition 



4.1 Defining the characteristic size 

To analyse the domain growth quantitatively in the following simulations we obtain the first zero crossing 
of the coordinate-space pair correlation function. This defines the characteristic domain size R{t). The pair 
correlation function is defined by: 

C(r,t)^^'^^"'^^^^^^ + '-'^)""' (41) 

where q(x, t) is the colour charge at site x, and the integral is taken over the whole system. We compute 
C(r, t) for each of our simulations, and perform an average over an ensemble of initial conditions. Taking 
the spherical average of C(r, t) yields C(r, i), the first zero of which gives the characteristic domain size. We 
obtain the first zero by performing a linear interpolation between the last point greater than zero and the 
first point less than zero. 



4.2 Scaling analysis 



The first test of scaling is applied to the correlation fmictions of Section 4.1. If the domain structures are 



self-similar at different times during the coarsening, the scaling function fix) defined as 



(42) 



should be independent of t. Alternatively the scaling criteria may be applied to the Fourier transform of 
C(r, i), the structure factor: 



5(k,t) 



1 



/ (7(x, t) exp (— ik • x)(ix'^ 



(43) 



Spherically averaging this quantity yields S(k,t), which has a similar scaling form: 



S{k,t)^L-^g{kL) = L-^g{q), 



(44) 



where g{q) is the fourier transform of /(a;). The functions f{x) and g{x) for simulations with /3 > 0.07 
are plotted as a function of x and q respectively in figs. || and ^. These datasets scale quite well, whereas 
datasets from simulations with (3 < 0.07 do not. The specific form of g(q) has been the subject of much 
theoretical attention. The three features which receive the most attention are the small, intermediate and 
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large q regimes. We show in figure |j the behaviour of our model in these three regimes. In the large q region 
we see a clear Porod tail with g{q) cx q^^ . This region ends when q probes the scale of the interface width. 
Furukawa speculated that there would be a regime g{q) oc q^^: we see a behaviour closer to g{q) oc (7^^, in 
agreement with Jury et. al. In the small-q limit we see g{q) oc q~^ ^ with 2 < (5 < 4. In ||2^ the constraint 
(5 < 4 was derived for a Cahn-Hilliard theory without hydrodynamics. However, this derivation assumed a 
dynamical scaling exponent of 1/3 and did not include fluctuations. Furukawa speculated that fluctuations 
would cause the g{q) oc q^ — g{q) cx q^ crossover at small q. 

The above analysis of the correlation functions and structure factors is independent of the form of R{t). 
We now introduce the reduced length and time variables I and t as defined in equations (|l^) and (^3|). To 
exclude finite size effects we use data for which the characteristic domain size is less than i of the system 
size. Data from our simulations satisfying this constraint spans a range of 3 < i < 152 and 1 < ^ < 10. A 
study using the DPD algorithm reached a range of 1 < Z < lO'^ and 1 < i < 5 x 10'*, and a latticc-Boltzmann 
algorithm has been used to reach a range of 1 < Z < 2 X 10^ and 1 < t < 5 X 10"^ ||l|,|l. The analysis method 
used to plot the data displayed in figures ^ and ^ from simulations with different values of Lq ^nd Tq is 
identical to that of Kendon et al. . 

We find that data with 0.10 < 13 < 100 show growth with an inertial (/ oc t^/'^) exponent, and R{t) 
collapses well onto a single scaling function. Simulations with [3 = 0.03 show slow growth where I oc t^/^ . 
Visualisation of the order parameter for this data shows that sharp interfaces do not form during the 
simulation, and so data for these low values of surface tension do not satisfy the criteria for the postulated 
universal scaling regime. This is confirmed by our scaling analysis of the correlation functions above, where 
data for 0.03 > (3 > 0.07 did not collapse onto the same curve shown in figure 3. The time evolution of the 
domain size in reduced units for 0.03 > (3 > 0.08 is shown in figure |^. 

This inertial exponent is not consistent with the previous results of Kendon et al. and Jury et al. 
although it is similar to the results of Appert et al. [||. In the (/ = ct^/"^) scaling of Appert et al. was 
ascribed to 'excessive diffusion' in the lattice gas algorithm. The particulate noise characteristic of the lattice 
gas is not present in their lattice-Boltzmann algorithm, and so it is perfectly plausible that these fiuctuations 
continue to inhibit phase separation, reducing the growth exponent in the viscous regime to some effective 
exponent close to the 2/3 usually associated with the inertial regime. An identical effect is well known in two 
dimensions, where models with fiuctuations yield an exponent of 1/2 in the early time regime, and models 
without yield an exponent of 1/3 p^,|4C|]. 

This issue deserves some closer examination, however. If we interpret our structure factor and corre- 
lation function data as exhibiting goood scaling collapse (i.e. for all times considered the morphology is 
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characterised by a single lengthscale) , then the g{q) oa — g{q) (x q^ crossover at small q is consistent with 
speculations by Furukawa that fluctuations would cause the q^ behaviour. As noted above these fluctuations 
may act to inhibit the phase separation and give rise to a lower effective exponent. 

Our immiscible lattice-gas model reduces to a simulation of a convection-diffusion equation for (3 = 0.0. 
The diffusion constant in this equation is a function only of the collision rules. It may be possible by a 
judicious choice of collision rules or addition of rest particles to vary this diffusion constant independent of 
surface tension to fully investigate its effect. 

If. however, we interpret the small q behaviour as being indicative of poor scaling collapse, the exponent 
seen may be interpreted as a crossover between an early time and viscous regime. This issue could only be 
resolved by larger and longer simulations, and present hardware limitations mean that this must remain a 
matter for future investigation. 

There is at present considerable uncertainty about both the universality of the asymptotic scaling and 
what the true asymptotic regime is. Theoretical work concerning the dispersion relation on fluid interfaces 
casts doubt on the universality of hydrodynamic phase separation in three dimensions p9[ . Recent work 
by Kendon | |39| ] proposes a new set of macroscopically deduced scaling laws, based on the energy balance 
in the system. This analysis suggests that there may be more than one length scale of importance in 
spinodal decomposition. Another derivation of the growth laws proposes that the n = | is transient and true 
asymptotic scaling is n = | p^ . However, it would be difficult to resolve such a subtle difference in exponent 
with any current numerical method. In addition, the derivation of p8[ | assumes a droplet morphology, even 
in the symmetric case. No such morphology is seen in any simulation of the symmetric case of which we 
are aware. The existence of multiple length scales and breakdown of scaling in two dimensional spinodal 
decomposition is well established [^,^, and the existence of a single underlying scaling function in three 
dimensions remains an open question. 



5 Kinetics of phase separation for off-critical binary fluids 

In this section we analyse the dynamics of domain growth in systems where the composition is asymmetrical. 
In such systems the domain structure is not bicontinuous. The minority oil (or water) phase exists as 
droplets, and the domains coarsen by diffusion of oil (or water) from the bulk onto the droplet, and by 
droplet coalescence. The composition (p for a system where oil is the minority phase is defined as: 

^=^^^ (45) 

Po + Pw 
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where po and are the reduced densities of oil and water respectively. We performed simulations for 
values of = 0.2, and — 0.4. These simulations were performed on 64"^ systems for 3000 time steps. 
The correlation function was averaged over five independent simulations for (p = 0.2, and three independent 
simulations for 4> ~ 0.4. The results are shown in figures ^ and ^. 

The data in our simulations gives effective exponents of 0.543±0.01 and 0.573±0.003 for compositions 0.2 
and 0.4 respectively. We do not observe the exponent 1/3 expected from simple theories of droplet coalescence 
and nucleation. However, as fig. |^ shows, the morphology for both these volume fractions is far from being 
a collection of isolated drops. It is likely that the exponent we measure is therefore intermediate between 
the critical case and droplet coalescence and nucleation. The existence of such intermediate exponents is 
well established. Previous work by Appert et al. [|j on off-critical decomposition for (j) ~ 0.05 also saw 
growth more rapid than the t^/^ expected for simple nucleation and coalescence. The work of for 
the two-dimensional implementation of our model similarly saw a continuous variation of exponent with 
composition, as did pij for equal viscosity DPD fluids. 



6 Self assembly in ternary amphiphilic fluids 

In this section we turn to the analysis of ternary amphiphilic fluids. We concentrate on systems with equal 
amounts of oil and water, varying the amount of surfactant for each simulation. The presence of surfactant 
dramatically alters the interfacial energetics and structure, and consequently the dynamics of domain growth. 
Specifically, the adsorption of surfactant at oil-water interfaces and its concomitant lowering of the surface 
tension weaken the forces which drive binary immiscible phase separation. For sufficiently large surfactant 
concentrations the final equilibrium state is a bicontinuous or sponge microemulsion phase, where the domain 
growth is arrested at some final characteristic domain size Re- Such a system is depicted in fig. (|l0[). All three 
components arc bicontinuously connected, with the surfactant particles lining the interface in a monolayer 
between the percolating oil and water regions. 

We work with (3 — 1.0 and use the values for the coupling coefficients in eqn (|2^) as defined earlier (in 
Sec. H) . The results that follow for the ternary emulsion system have been obtained from a 64'^ system with 
periodic boundary conditions in all directions, the system having been intialised from a random quench with 
the particles placed randomly on the lattice. The total reduced density of the system was kept constant at 
0.5. The characteristic domain size R{t) was measured from the spherically averaged pair correlation function 



as described in Section 4.1. The growth of sponge as opposed to droplet phases in these systems means that 



for the lowest values of the surfactant concentration the growth behaviour should represent a perturbation of 
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the critical quenches investigated in Section ^. The results for reduced surfactant concentrations 0.02, 0.06, 
and 0.08 support this conclusion. In all cases the late-time behaviour is consistent with an algebraic growth 
law of the form R ^ t^l'^ . Prior to this there is an an early-time regime in which the growth is consistent 
with a diffusive algebraic exponent of i. Visualisation of the surfactant densities during this period shows 
that surfactant adsorbs at the oil-water interfaces. 

Once the system reaches the biconinuous oil and water state the inertial hydrodynamic regime begins. 
The length of the diffusive period increases with increasing surfactant concentration (50 time steps for 
reduced surfactant concentration = 0.02, 60 time steps for ps = 0.04, 70 time steps for p^ = 0.06 and 100 
time steps for ps — 0.08). As one increases the surfactant concentration the fluctuations in the domain size 
at late times increases. This is consistent with the known dynamical nature of the adsorption and desorption 
of surfactant to and from interfaces in ternary systems as one approaches emulsification. With a reduced 
surfactant density of 0.10 the time evolution of the domain size ceases to follow the algebraic R ~ t^/'^ 
growth law. We see a slowcr-than-algebraic growth law, as previously observed in the two dimensional 
implementation of our model where a growth fimction R{t) = (Inf)^ was proposed, based on a comparison 
with slow growth in systems with quenched impurities ||3^, Consequently we look at a plot of hit 
against domain size in order to determine whether we have logarithmically slow growth in this regime. The 
characteristic domain size R{t) for surfactant concentrations 0.10,0.12,0.14,0.16 are plotted against hit on 
logarithmic scales in fig. (p^). One can clearly see a transition through a regime where logarithmically slow 
growth dominates throughout the timcscalc of the simulation. Such a growth law is inconsistent with arrest 
of the domain growth at late times for reduced surfactant density below 0.16. As we increase the reduced 
surfactant density beyond 0.16 we see complete arrest of the domain growth at late times. We performed a 
fit of a stretched exponential function to datasets from simulations with ps = 0.18. This function is defined 
identically with : 

R{t) = A-B exp(-Ct^) (46) 

These fits are shown in fig. dl^). To quantify the effect that the surfactant has on the domain size at late 
times in these simulations, following Emerton et al. and Gyerc et al. [i^ , we define Rc as the domain 
size at time step 850, where data which is unaffected by finite size effects is available for all surfactant 
concentrations. We plot Rc against the inverse of the reduced surfactant density at the interface (having 
subtracted the miccUar and monomcric concentration) in fig. ^ (for systems with quenched impurities the 
amplitude of the disorder is analogous to the surfactant concentration). We find a linear dependence of Rc 
on l/ps, consistent with the results found in [Q and |^]. 



17 



7 Conclusions 



We have studied the dynamics of domain growth in binary immiscible fluids, in both critical and off-critical 
cases, and in ternary (oil-water-surfactant) emulsions and microemulsions, using a three-dimensional hydro- 
dynamic lattice gas model. In the critical binary case we performed a scaling analysis similar to that of 
Jury et al. ^ and Kcndon et al. and find late-time growth behaviour which is consistent with an inertial 
hydrodynamic exponent. However the position of the late-time domain growth on a plot of reduced time 
variables is not consistent with that of in the same way that the results of Appert et al. 

showed an inertial exponent in a crossover region of the reduced scaling plot. In the off-critical case we find 
effective exponents of 0.57 and 0.54 for compositions (j) = 0.4 and 0.2 respectively. In the microemulsion 
case we observe a slowed algebraic growth with exponent |, followed by a regime in which we see evidence 
for logarithmically slow growth. For concentrations of surfactant high enough to completely arrest domain 
growth we observe good agreement with a stretched exponential fit to our data. Overall, all our results are 
fully consistent with behaviour we observed previously in our two-dimensional lattice gas model. 
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Figure 1: Surface tension as a function of inverse temperature-like parameter (3 (both in lattice units) for a 
binary immiscible fluid. Error bars are not shown, but arc smaller than the symbols. System sizes are 64'^. 
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Figure 2: Surface tension (lattice units) as a function of surfactant concentration in a ternary system. Squares 
are values calculated after equilibration, triangles are calculated as a function of instantaneous surfactant 
density. All simulations have /5 = 1.0 and were run on 64 x 32 x 32 lattices. 
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Figure 3: The scaling function f{x) as defined in eqn (|4^) from the pair correlation function in critical 
binary phase separation for 0.08 < /3 < 100 (lattice units). Data is taken from all simulations in this range 
of inverse temperature-like parameter /?. System sizes are 64'^ and 128'^; data is taken between 100 and 2300 
time steps. 
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Figure 4: The scaling function g{x) as defined in eqn (Q) from the structure factor in critical binary phase 
separation for 0.08 < f3 < 100 (lattice units). Data is taken from all simulations in this range of inverse 
temperature-like parameter (3. System sizes are 64^ and 128'^; data is taken between 100 and 2300 time 
steps. The dotted line is {kL?'), the long dashed line is {kL'^), the short dashed line is {kL~'^) and the solid 
line is {kL~'^). They are included as guides to the eye. 
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Figure 5: Scaling plot in reduced variables {l,t) for critical binary immiscible fluid phase separation. Data 
are from simulations with /3 = 0.03, 0.04, 0.06, 0.08 (from left to right). System sizes are 64'^. 
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Figure 6: Scaling plot in reduced variables (l,t) for critical binary immiscible fluid phase separation. Data 
are from all simulations with 0.10 < (3 < 100. System sizes are 64^ and 128^; data is taken between 100 and 
2300 time steps. The solid line has gradient 2/3 and is included as a guide to the eye only. 
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Figure 7: Interface morphology at time step 500 for minority phase (oil) volume fraction 0.20 (a) and 0.40 
(b). System size is 64'^. 
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Figure 8: Scaling in off-critical binary phase separation for /? = 1 with composition </) = 0.2 (lattice units). 
The solid line is a least squares fit with effective exponent 0.54±0.01. Error bars show one standard deviation 
on the mean of five independent simulations. 
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Figure 9: Scaling in ofF-critical binary phase separation for = 1, at a composition 4> ~ 0.4 (lattice units). 
The solid line is a least squares fit with effective exponent 0.573 ± 0.003. Error bars show one standard 
deviation on the mean of three independent simulations. 




(a) (b) (c) 

Figure 10: Sponge microcmulsion phase at time step 850 following a random initialisation, (a) Water slice 
plane, (b) Surfactant slice plane, (c) Oil slice plane. The system size is 64'^. Reduced densities of oil, water, 
and surfactant in this system are 0.19,0.19 and 0.12, respectively. 



28 



1 I , , , I 

100 1000 

Tine (time steps) 

Figure 11: Time evolution of characteristic domain size (lattice units) for a ternary amphiphilic fluid with 
reduced densities of oil, water and surfactant 0.24, 0.24 and 0.02 respectively. The solid line is a fit with 
an effective exponent of 0.59. The error bars show one standard deviation on the mean of five independent 
simulations. The system size is 64'^. 
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Figure 12: Crossover in ternary domain growth behaviour to logarithmicaUy slow growth, and then to an 
intermediate regime. Surfactant concentration 0.10 shows behaviour in a regime crossing over from algebraic 
to logarithmic growth. Surfactant concentrations 0.12 and 0.14 show convincingly logarithmic growth, while 
concentration 0.16 shows time evolution in a regime between logarithmic and stretched exponential growth. 
Both time and characteristic size are measured in lattice units. 
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Figure 13: Plot of characteristic domain size against time (both in lattice units) for a ternary system with 
reduced surfactant density 0.18. The solid line is a least squares fit of a stretched exponential function to 
the data. Error bars show one standard deviation on the mean of five independent simulations. The system 
size is 64^. 
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Figure 14: Plot of final characteristic domain size Rc (lattice units) against the inverse of the reduced 
surfactant density of surfactant 1/ps in the system. We have corrected ps by removing the micellar and 
monomeric surfactant concentrations. 
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